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The identification of trajectories that contribute to the reaction rate is the crucial dynamical ingredient in any classical 
chemical reactivity calculation. This problem often requires a full scale numerical simulation of the dynamics, in 
particular if the reactive system is exposed to the influence of a heat bath. As an efficient alternative, we propose here to 
compute invariant surfaces in the phase space of the reactive system that separate reactive from nonreactive trajectories. 
The location of these invariant manifolds depends both on time and on the realization of the driving force exerted by 
the bath. These manifolds allow the identification of reactive trajectories simply from their initial conditions, without 
the need of any further simulation. In this paper, we show how these invariant manifolds can be calculated, and used in 
a formally exact reaction rate calculation based on perturbation theory for any multidimensional potential coupled to a 
noisy environment. 
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I. INTRODUCTION 

Transition State Theory (TST) provides the conceptual 
framework for large parts of reaction rate theory. Original ly 
developed to describe the reactivity of small molecules*^ it 
was later extended to encompass a wide variety of processes 
in different branches of science, whose only commonality is a 
transition from well-defined "reactant" to "product" states^— 
The reason for this success is that TST proposes a simple 
answer to the two central problems of reaction dynamics: It 
identifies a reaction mechanism, and provides at the same time 
a simple approximation to the reaction rate. 

More specifically, TST is based on the observation that the 
rate limiting step in many reactions is the crossing of an en- 
ergetic barrier. The top of this barrier then forms a bottleneck 
in the phase space of the reactive system. A reaction can only 
take place if the barrier is crossed. If a dividing surface (DS) 
between reactant and product regions of phase space is placed 
close to the bottleneck, the reaction rate can be computed from 
the steady-state flux through that surface. A strictly recross- 
ing free DS can be constructed in the phase space of reac- 
tive systems with arbitrarily many degrees of freedo m 12 ' 14,15 . 
The simplest approximation to the rate is then obtained under 
the assumption that reactive classical trajectories cross the DS 
only once and never return. This assumption is often appro- 
priate for reactions in the gas phase if the DS is adequately 
chosen, but even then many reactions strongly violate this as- 
sumption. Moreover, if the system is strongly coupled to an 
environment, for example a liquid solvent, the no-recrossing 
assumption is usually impossible to enforce strictly, and often 
any DS is crossed many times by a typical trajectory. As a 
result, a TST rate calculation significantly overestimates the 
reaction rate. For this reason, the focus of TST has long been 
to construct a DS that eliminates or at least minimizes recross- 
ings (see Ref.[l6|for a review). 

The recrossing problem can be solved if the reactive trajec- 



tories that contribute to the rate can be identified reliably. An 
obvious means to this end is the numerical simulation of rep- 
resentative trajectories under the influence of the environment. 
However, such calculations are usually very time consuming . 
The advantage of the TST approximation is its simplicity. It 
identifies reactive trajectories simply by noting that they cross 
the DS from the reactant to the product side. This criterion, 
which fails if recrossing cannot be ruled out, is easy to use be- 
cause it only takes account of the instantaneous velocity with 
which a trajectory crosses the DS. Nevertheless, it raises the 
prospect of a criterion to identify reactive trajectories simply 
from their initial conditions, without the need to study their 
time evolution. In the present paper we will derive such a cri- 
terion and demonstrate how it can be used in a rate calculation. 

The Langevin equation has been widely used to model 
the interaction of a reactive system with a surrounding heat 
bathfii^ Being a classical model, this description neglects 
quantum effects such as barrier tunnelling, which can be im- 
portant in the case of light particles 2 ^, and the interaction 
with excited surfaces through conical intersections 21 . In this 
setting, Kramers* 2 - explicit derived expressions for the rate 
of escape across a parabolic barrier that apply separately in 
the limits of weak and strong damping. The generalized 
Langevin equation is equivalent to a Hamiltonian model in 
which the reactive system is bilinearly coupled to a bath of 
harmonic oscillators^ This reformulation allowed extensions 
of Kramers' rate theory that apply to situations with arbi- 
trary frictio n 24,25 or that include corrections due to anhar- 
monic barriers^— In this respect, it has long been predicted 
that the rates of activated processes should rise with the cou- 
pling to the solvent in the weak coupling regime. However, 
its direct observation in particle-based models had been elu- 
sive because the coupling typically places the processes in the 
spatial-diffusion limited regime wherein rates decrease with 
increasing friction. Recently, the Kramers turnover in the rate 
with microscopic friction has been observed in molecular dy- 
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namics trajectories calculation of the LiNC^LiCN in a bath 
of Ar atoms^2 This observation provided direct and unam- 
biguous evidence for the energy-diffusion regime in which 
rates increase with friction. In the present work we will not 
consider any explicit Hamiltonian model for the heat bath; its 
influence will instead be described by means of a Langevin 
equation. This approach allows to work within the finite- 
dimensional phase space of the reactive system alone, rather 
than the infinite-dimensional phase space of the bath. This 
is advantageous from a computational point of view and also 
conceptually convenient because the phase space is easier to 
visualize in low dimension. 

The aim of this paper is to describe the geometric phase 
space structures that allow to classify a trajectory as reac- 
tive or nonreactive just by looking at its initial condition, thus 
avoiding the need of carrying out a numerical simulation. Be- 
cause the fate of a trajectory with a given initial condition de- 
pends on the external force to which it is exposed, any such 
criterion must take account of the precise realization of that 
force. A general framework to do that was proposed in a re- 
cent series of papers?22~~ 4 including the identification of re- 
active trajectories^ and the rate calculation^. It was there 
shown that the Langevin equation gives rise to a specific tra- 
jectory called the Transition State (TS) trajectory that remains 
in the vicinity of the energetic barrier for all times, without 
ever descending into any of the potential wells. This TS trajec- 
tory depends on the realization of the noise, and takes over the 
role of the fixed saddle point in the conventional TST A cru- 
cial observation in Refs. [3(3 and HI] for the case of a harmonic 
barrier is that the dynamics described by the Langevin equa- 
tion become noiseless when expressed in a time-dependent 
coordinate system for which the TS trajectory is the moving 
origin. In the system of relative coordinates it is easy to iden- 
tify a TST DS that is rigorously free from recrossing. It gives 
rise to a DS in the original, space fixed coordinate system 
that is still recrossing-free. This DS is time-dependent since 
it is attached to the TS trajectory, and it moves through phase 
space with it. Even more significantly, this construction yields 
surfaces in phase space that separate reactive from nonreactive 
trajectories. These surfaces are the stable and unstable man- 
ifolds of the TS trajectory, and they also depend on time and 
on the realization of the noise. Once they are known, initial 
conditions on one side of the surface are immediately classi- 
fied as reactive, while those on the other side are nonreactive. 
Thus, the existence of these invariant manifolds solves the di- 
agnostic problem of standard rate theory that was explained 
above. They were used in Ref. [33] to obtain a compact rate 
formula, strictly valid only for harmonic barriers. An ad hoc 
application to systems with an anharmonic barrier produced, 
however, promising results^^ 

In the present paper, we develop a rigorous generaliza- 
tion of the time-dependent TST formalism applicable to an- 
harmonic barriers using perturbation theory. We show that 
the invariant manifolds persist in anharmonic systems and, 
more importantly, they retain the ability to distinguish be- 
tween reactive and nonreactive trajectories, thus determining 
the chemical reactivity of the system. Finally, a simple pertur- 
bative scheme that allows one to calculate the invariant man- 



ifolds for a specific anharmonic potential barrier will be pre- 
sented, and it will be used to obtain an analytic expansion for 
the reaction rate. In the first part of the paper, we restrict our 
study to the one-dimensional case. In this situation, the finite 
barrier corrections that were obtained in Refs. I261 - I281 will be 
recovered. We have already given a brief account of these re- 
sults in Ref. EL We will here supply the details of the calcula- 
tion that could not be presented within the confines of a Com- 
munication. We will then introduce the modifications to the 
theory that are necessary to accommodate multidimensional 
reactive systems. The efficacy of our method is demonstrated 
by deriving the first and second order corrections to the reac- 
tion rate in the two-dimensional model potential already used 
in Refs. US and[M 



A final point is worth commenting on in this Introduction. 
Perturbative rate calculations on multidimensional anhar- 
monic barriers have also been recently reported in Refs. l36l - 
l39i As in the present work, these authors based their work on 
the identification of the TS trajectory for the harmonic limit in 
Refs.^anddll Our work, however, goes beyond those previ- 
ous results in two main respects. First, and most importantly, 
it provides an explicit and detailed description of the invari- 
ant geometric structures in phase space that govern the reac- 
tion dynamics, rather than studying them implicitly through 
approximate invariants and their imprint on an ensemble of 
trajectories. Second, whereas the normal form procedure in 
Refs. 13614391 aims at constructing a coordinate system in which 
the dynamics in the neighborhood of the barrier can be simpli- 
fied in general terms, we derive a version of the perturbation 
theory that is specifically directed at calculating the invariant 
manifolds that are relevant to reaction rate theory. This pertur- 
bative scheme can therefore be much simpler, and permits the 
analytical computation of corrections to Kramers' transmis- 
sion factor for anharmonic potentials. Indeed, the calculation 
of the invariant manifolds can be easily carried out by hand, 
whereas a normal form transformation always requires com- 
puter assitance. This ease of computation makes the invariant 
manifolds an attractive tool for practical rate calculations. 



The outline of the paper is as follows. In Section [TT] we 
present the basic definitions and results of rate theory that 
will be used to develop our method. Section [III] is devoted 
to a qualitative description of the invariant manifolds that give 
structure to the dynamics in the vicinity of an energy barrier, 
and section|IV]presents a method for their calculation. In sec- 
tion [V] a general expression for the reaction rate in the case 
of an anharmonic barrier is derived. A description of the sta- 
tistical properties of the invariant manifolds that are required 
to evaluate the rate formula, the perturbative and numerical 
results for various one-dimensional potentials are also given. 
Finally, in section[VT]we discuss the modifications to the fore- 
going developments that are required in multidimensional sys- 
tems, and we also present results for the reaction rate on an 
anharmonic two-dimensional barrier. 
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II. FUNDAMENTALS OF RATE THEORY 

In this section we outline the fundamentals of reaction rate 
theory that will be used in the rest of the paper. The reader is 
referred to Ref s . [l7l - fT9l for more details. 

We assume that the reactant and product regions in configu- 
ration space are separated by a DS that is characterized by the 
value x — x* of a generalized reaction coordinate x, which we 
choose such that the product region is given by x > x*. The 
reaction rate is then given by the flux-over-population expres- 
sion 



N 



(1) 



where N is the average population of the reactant region and 



J = (VxXi(y x ,q x ,Vi.)) a i C 



(2) 



is the reactive flux out of that region. Here, v x denotes the ve- 
locity component perpendicular to the DS, q ± the coordinates 
within the surface and v ± the corresponding velocities. The 
characteristic function Xr( v x>q±>V±) takes the value 1 if the 
trajectory starting at x — x*,v x ,q ± ,v± is reactive, i.e., moves 
to products for large times, and otherwise. Its purpose is to 
ensure that only reactive trajectories contribute to the reactive 
flux. The average in Eq. (f2]i extends over the realizations, a, 
of the external noise and over a thermal equilibrium ensemble 
of initial conditions that are constrained to lie on the DS. The 
latter ensemble is described by a probability density function 



p(x, v x , q ± , v x ) = 5{x - x*) exp P^in v ± 



), 0) 



which includes a Boltzmann distribution of the velocities v x 
and a Boltzmann distribution 

1 / v 2 J2 + U{xK qi ) \ 
Px(qx,vx) = r xv[ — j (4) 

of the transverse coordinates and velocities. The factor Z in 
Eq. (|4]i is the partition function of the transverse motion. It 
ensures that 



dq ± dv ± p(q ± ,v ± ) = 1- 



In Eq. ([3} we have used mass-scaled coordinates and we have 
left out an overall normalization factor. In particular, we did 
not include the Arrhenius factor 



exp 



that includes the activation energy AE* of the reaction. The 
overall normalization of the distribution function is well un- 
derstood, and it is irrelevant to the calculation of the transmis- 
sion factor ^ below, on which we will focus in this work. 
For simplicity, we can therefore work with the unnormalized 
distribution function 



The characteristic function % t in Eq. (O encodes the entire 
complexity of the reaction dynamics on an anharmonic bar- 
rier. The main task of a reaction rate calculation is to evaluate 
this function. In general, this can only be achieved by a nu- 
merical simulation. A simple approximation to this crucial 
ingredient is provided by TST. It assumes that no trajectory 
can cross the DS more than once. As a consequence, every 
trajectory that crosses the DS from the reactant to the prod- 
uct side must be reactive, every trajectory that crosses in the 
opposite direction must be nonreactive. To implement this ap- 
proximation, we replace the characteristic function in Eq. (f2| 
by 



X TST (v x ,q x ,v x ) = 



Vx > 0, 
v x < 0. 



(5) 



This gives rise to the TST approximation to the rate constant 
jTST {v xX TST (v x ,qx,vj) 



lIC 



N 



(6) 



in which the average over the noise a can be suppressed be- 
cause x TST does not depend on it. 

When the no-recrossing assumption of TST is not satisfied, 
the approximation (0 will overestimate the rate, often by a 
large factor. To quantify the effects of non-TST behavior, a 
transmission factor, 

k 

K - < 1 

£TST - ' 

is introduced that relates the exact rate to the TST approxima- 
tion. It can be obtained from the ratio of the flux across the 
barrier to its TST approximation: 

(v^l(v. t ,g_L,V_L)) or ,IC (?) 

(v x x TST (vx,q±,v ± )) lc 

To evaluate © numerically, one can randomly sample initial 
conditions and noise sequences from the appropriate ensem- 
bles, and simulate the behavior of each trajectory until its en- 
ergy is so far below the barrier top that it can be regarded as 
having been thermalized on either the reactant or the product 
side of the barrier. The trajectory can then be classified as re- 
active or non-reactive depending on what state it reached. All 
numerical results presented in this work were obtained in this 
way. 

This algorithm is conceptually straight-forward, but com- 
putationally costly. It would be highly desirable to find a cri- 
terion that allows one to identify the reactive trajectories with- 
out having to carry out a numerical simulation. The following 
sections will describe the phase space structures that will pro- 
vide such a criterion. 



III. TIME-DEPENDENT INVARIANT MANIFOLDS 

A. The Langevin model 

We begin by specifying the model that will be used. The 
Langevin equation describes the reduced dynamics of a low- 
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dimensional system coupled to an external heat bath^ It is 
given by 



-V q U(q)-Tq+£ a (t), 



(8) 



where q is an N dimensional vector of mass-scaled coordi- 
nates, U(q) is the potential of mean force, T is a symmet- 
ric positive-definite N x N matrix of damping constants, and 
€ a (t) is the fluctuating force exerted by the heat bath. It is con- 
nected to the friction matrix T by the fluctuation-dissipation 
theorem^ 



(f a (t)? a (t')} a = 2k B TT6(t-t'), 



(9) 



where ks is the Boltzmann constant and T is the temperature. 
Throughout most of this work, we consider a one-dimensional 
problem in which the friction matrix Y simply reduces to a 
scalar y, and the position vector q contains a single coordinate 
x. If we expand the potential of mean force around its saddle 
point, we can write it as 

U(x) = -\oj\x 1 + eyx 3 + s 2 ^/ + . . . . (10) 

where s is a formal perturbation parameter that serves only to 
keep track of the orders of perturbation theory, and finally will 
be set to s — 1 . For the mean force itself we write 



dU 
dx 



u b x + f(x), 



(11) 



where fix) denotes the anharmonic parts of the force. 



B. Time-dependent transition states 

Because the Langevin expression d8) is a second order dif- 
ferential equation, its phase space is two-dimensional, with 
coordinates x and v x = x. As it was observed in Refs. 
[32l the dynamics of the Langevin equation in the harmonic 
approximation can be diagonalized by rewritting it in coordi- 
nates u and s given by 



u — 



v x ~ ^ s x 
A u - /l s 



~Vx + *v.X ,,,,, 

s=— — , (12) 



or 



x — u + s, 
The constants 



v x = A u u + AsS. (13) 



A s , u = -i(r± ^r 2 + 4^) (14) 



are the eigenvalues that arise in the diagonalization. They sat- 
isfy A s < < A a and 

A u + As = -y, A U A S = -a%. 

In the new set of coordinates, the equations of motion read 

fix) 1 

U = A U U H 1 %ait), 

Aa - A s A u - A s 



s — ArS — 



1 



fix) 

A u — A s A a — A s 



?o(t). 



(15) 



These equations decouple in the harmonic approximation, i.e., 
if fix) = 0, but they are still subject to the time-dependent 
stochastic driving force This time dependence can be 

removed by the coordinate shift 



Au — u — u 



As 



(16) 



where 

k*(0 



-!— S[^,?„; t], s*it) = -— !— S [As,£ a ; f], 
— A« A u — A s 

(17) 



and the S functional a 30 - 41 are given by 

gir) exp(//(f - t)) dr : Re /i > 0, 



Sr\M,g;t] 



V \J —CO 



(18) 



(t) exp(ju(f - t)) dr : Re fi < 0. 



The subscript t is used in the S functional to indicate the inte- 
gration variable. This subscript will be left out whenever this 
does not cause any ambiguities. Similarly, we have for the 
sake of simplicity not indicated in our notation that M*(f) and 
s*(f) depend on the realization a of the noise, although they 
both obviously do. 

The functions M T (f) and solve the equations of mo- 
tion in the harmonic limit fix) = 0. They can therefore be 
regarded as the coordinates of a special trajectory called the 
TS trajectory. This trajectory is distinguished from all other 
trajectories that are exposed to the same noise by the fact that 
it remains in the vicinity of the saddle point for all times, 
whereas a typical trajectory would descend into either the re- 
actant or the product well both in the remote past and in the 
distant future. Accordingly, when using coordinates Au and 
As, we are describing a trajectory relative to the TS trajectory, 
which acts as a moving coordinate origin. In what follows, 
we will refer to Au and As as relative coordinates and to the 
original u and s, or x and v x as space fixed coordinates. 

The equations of motion in relative coordinates are 



Aii — A u Au + 
As = A s As - 



fix) 
fix) 

A u — /Is 



(19a) 
(19b) 



At first sight, it appears that the time-dependent and stochas- 
tic shift ( TTol ) has removed both the time-dependence and the 
dependence of the realization a of the noise. However, this 
is only true in the harmonic approximation. If we express the 
position coordinate x in terms of the relative coordinates Au 
and As as 



x = x* + Au + As, 
with a* = k* + s*, Eq. ( fT9l ) turns into 

fix* + Au + As) 



Aii = A u Au + 



A u — /is 



(20) 



(21a) 
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f(xr + Au + As) 
As = AsAs-— -. (21b) 

The position x*(t) of the TS trajectory represents a time- 
dependent stochastic driving in these equations of motion. 
Nevertheless, the coordinate shift has removed the stochas- 
tic driving from the leading-order terms in (OTl and pushed it 
into the anharmonic perturbation. 

The description of the geometric phase space structure in 
the vicinity of the saddle point is most easily done if one starts 
from the harmonic limit. A full discussion can be found in 
Refs.l3(il and [3ll The equations of motion jl9[ decouple and 
become time independent when f(x) = 0, and they can then 
be easily solved by writing 

Au(t) = Am(0) e A "', 

As(t) = As(0)e**'. (22) 

Since A u > and A s < 0, the coordinate Au grows expo- 
nentially in time, whereas As shrinks. Therefore, Au and As 
correspond to unstable and stable directions in phase space, 
respectively. In particular, the lines Au = and As = are in- 
variant under the dynamics. A trajectory that starts on the line 
Au — will asymptotically approach the origin as t — > 00; this 
line is called the stable manifold of the origin. A trajectory on 
the line As = will move away from the origin as t — > 00, but 
it will approach the origin as t — > -00; this line is called the 
unstable manifold of the origin. 

The stable and unstable manifolds of the origin, together 
with several typical trajectories in relative coordinates, are 
shown in Fig. QJ. The invariant manifolds separate trajec- 
tories with different qualitative behavior. Trajectories above 
the stable manifold, i.e., with larger relative velocity, move to 
the product side of the barrier for asymptotically long times, 
whereas trajectories below the stable manifold move to the 
reactant side. Similarly, trajectories above the unstable man- 
ifold come from the reactant side in the distant past, whereas 
trajectories below the unstable manifold come from the prod- 
uct side. 

For a reaction rate calculation we need to ascertain whether 
a trajectory will turn into reactants or products in the future. 
In our approach this sentence is rephrased into the condition: 
We need to decide whether a trajectory lies above or below the 
stable manifold. In other words, the stable manifold encodes 
the information about the reaction dynamics that is most rele- 
vant to us. We will therefore focus on the stable manifold in 
what follows, largely ignoring the unstable manifold. 

We can return to space fixed coordinates by undoing the 
time dependent shift ( TToT l. After the shift, the stable and unsta- 
ble manifolds are not attached to the origin of the coordinate 
system any more, but instead to the TS trajectory as a mov- 
ing origin, as shown in Fig. 02b). Since the TS trajectory is 
time dependent, the manifolds will move through phase space 
with it. Nevertheless, they still separate trajectories with dif- 
ferent asymptotic behaviors. Given a trajectory with a given 
initial condition at a certain time, it can be classified as reac- 
tive or non-reactive by knowing the instantaneous position of 
the stable manifold at that time. Through the TS trajectory, 




FIG. 1. Phase space view of the time-dependent invariant mani- 
folds of the Langevin equation, (a) Invariant manifolds are time- 
independent in the harmonic approximation and in relative coordi- 
nates, (b) In space-fixed coordinates, the invariant manifolds are at- 
tached to the TS trajectory and move through phase space with it. 
(c) Anharmonic coupling deforms the manifolds. Both their position 
and their shape are stochastically time dependent, (d) Invariant man- 
ifolds can deviate strongly from the harmonic approximation if the 
anharmonicities are strong. 

that instantaneous position will depend on the realization of 
the noise. 

It is clear from Fig. |TJb) that at any time and for any re- 
alization of the noise the stable manifold intersects the axis 
x — at a point with a velocity V*. Trajectories with initial 
positions x = and initial velocities v x > V* are reactive, 
while trajectories with initial velocities v x < V* are not. The 
critical velocity V* depends on time and on the realization of 
the noise. For the harmonic approximation, it was shown in 
Ref.HH and it will be rederived below, that 

V* = V* = (A - A s )uH0). (23) 

Since the critical velocity characterizes reactive trajectories, 
the transmission factor (jTj can be expressed in terms of V* 
(see Ref . [33l and Section[V]below). 

This picture of the invariant manifolds was introduced in 
Ref s. [30I and[3j] and applied to rate calculations in Refs. [32| 
and L3J. The main purpose of the present work is to explore 
how this picture changes when anharmonicities of the barrier 
potential are taken into account. In this case the equations of 
motion ( fT9b are coupled in a nonlinear time-dependent way, 
and they cannot be solved easily. However, as long as the 
coupling is sufficiently weak, it can be expected to find a TS 
trajectory and with its associated stable and unstable mani- 
folds that are close to those in the harmonic approximation. 
Indeed, there are general theorems in the theory of stochastic 
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dynamical systems that guarantee the persistence of these 
structures. As shown in Fig. Q3 C X the invariant manifolds in 
an anharmonic system will be tangent to their harmonic ap- 
proximations at the TS trajectory, but they will not be straight 
lines anymore. Because the coupling term in jl9\ is stochasti- 
cally time dependent, the shapes of the invariant manifolds as 
well as their positions in phase space depend on time and on 
the realization of the noise. 

The intersection of the stable manifold with the axis x — 
will give rise to a critical velocity V* such that trajectories 
with initial velocities larger than V* will be reactive, those 
with smaller initial velocities will not. The critical velocity 
can therefore be used in a rate calculation in an anharmonic 
system just as it can in the harmonic approximation, though 
its value will be different from (|23V A method to calculate the 
critical velocity will be developed in Section lTVl 

In general it cannot be guaranteed that there will only be a 
single intersection between the stable manifold and the axis 
x = 0. In fact, if the reaction potential has wells on the re- 
actant and/or product side of the barrier, it is likely that there 
will be further intersections, as illustrated in Fig. H}d). If a tra- 
jectory on the stable manifold is followed backwards in time, 
it will descend from the barrier, settling in one of the wells for 
some time. If it is followed for long enough, it will eventually 
cross the barrier again into the other well. In doing so, it must 
cross the line x = again, and thus give rise to additional 
intersections between the stable manifold and that line. How- 
ever, as these additional intersections stem from previous bar- 
rier crossings, they must be neglected in the rate calculation. 
Only for extremely strong nonlinearities additional intersec- 
tions that are not separated by periods in which the trajectory 
was equilibrated in one of the wells will be found. We will ne- 
glect that possibility in what follows. Instead, we will apply 
perturbation theory to calculate a value for the critical velocity 
that reduces to its harmonic approximation in the appropriate 
limit. 

The TS trajectory (fTT) solves the equations of motion JT5t 
in the harmonic limit, but not in the presence of anharmonic 
coupling. Strictly speaking, therefore, Eq. ([PTt does not de- 
fine a TS trajectory on an anharmonic potential. Such a tra- 
jectory could be obtained by a perturbative expansion similar 
to the one to be developed in Section [IV] For our purposes, 
however, this will not be necessary. The harmonic TS trajec- 
tory forms a suitable basis for the perturbation theory. We will 
therefore use the notation m*, s* and jc* exclusively to denote 
the harmonic approximation to the TS trajectory. 



IV. PERTURBATIVE CALCULATION OF THE STABLE 
MANIFOLD 

The critical velocity is defined by the intersection of the 
line x = with the stable manifold of the TS trajectory. The 
stable manifold contains all those trajectories that approach 
the TS trajectory as t — > oo. They remain bounded for large 
times. Solutions to the equations of motion ( I2TT ) that satisfy 
this boundary condition at large time lie on the stable mani- 
fold. 



Equation ( 12 lab can be formally solved in terms of the 
S functional ( fT~8T > as 

Au(t) = Ce' 1 "' + — - — S [A u ,f(x* + Au + As); t]. 

A u - A s 

Notice that this is only a formal solution due to the presence 
of the unknown function Am in the r.h.s. of the equation. Fur- 
thermore, the S functional is undefined for most trajectories, 
only existing for the trajectories that remain bounded in the 
remote future. However, these are precisely the trajectories 
we are interested in. For consistency, we must then set C — 0, 
just as was done in Refs.l30landl3l1in the construction of the 
TS trajectory. A trajectory on the unstable manifold therefore 
satisfies the integral equation 

Au(t) = — - — S [A u , f(x* + Am + As); t]. (24) 

A u - 

This expression automatically incorporates the boundary con- 
dition at t — > oo that we wish to impose. 

For the stable component, we might be tempted to use the 
analogous formal solution 

As(t) = Ce u — S v + Am + As); t]. 

A u - A% 

However, the S functional for a negative eigenvalue depends 
on the infinite past of its argument and is well defined only for 
trajectories that remain bounded in the past. Most trajectories 
on the stable manifold, except for the TS trajectory itself, will 
not satisfy this condition. This difficulty can be circumvented 
by using the modified S functional 

S r [M,g;t]= f g (T)e^'- T) dT (25) 
Jo 

that is well defined for all values of yu. It satisfies the differen- 
tial equation 

d - 

-S[M,g;t]= f xS[M,g;t]+g(.t) 
at 

and the initial condition S [p., g; 0] = 0. With this functional, a 
formal solution to the equation of motion ( 121bt can be written 

as 

As(t) = As(0)e A >' — S [A\, fix* + Am + As); f]. (26) 

A u — A s 

Note that this integral equation does not impose any boundary 
condition on the function As, thus leaving free choice of the 
initial condition As(0). 

The critical velocity V* is determined by the condition that 
the trajectory with initial conditions x{0) = and v(0) = 
satisfies the integral equations (l24l and d26l i. The first one of 
these conditions can be rewritten as 

As(0) = -jc*(0) - Am(0), 

such that the initial condition for As, which is needed in 
Eq. (|26| |. is known once the initial value of Am has been de- 
termined from Eq. (l24l . The critical velocity can then be ob- 
tained from 

y* = v (0) = A u u(0) + A s s(0) 
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= (A u - A s )u(Q) as x(0) = u(0) + s(Q) = 

= (A a -A s )[u*(0) + Au(0)]. (27) 

In the harmonic approximation the trajectory that starts in 
the DS x = and lies in the stable manifold is given by 

Ako(0 = and As (f) = -x*(0)e' lsf . (28) 

For this case, Eq. (|27| | leads back to the result i 

V* = (A - A s )M*(0). 



When the solution (1281 1 is substituted into the integral equa- 
tions ( f24b or d26l i. the coordinate x = + Aw + As is replaced 
by 



X(t) = xHt) - e A ''x*(0) 



(29) 



This function represents the harmonic approximation to the 
coordinate x{f) of the trajectory under study. Moreover, it con- 
stitutes a suitable basis of the perturbative expansion. 

The leading-order correction to the critical velocity can be 
obtained from d24t as 



At/lead(0 = 



1 



S[A u ,f(X);t], 



A u — A s 
from which it follows that 

V* = S[A a ,f(X);0] 



(30) 



To obtain higher-order corrections to the critical velocity in 
a systematic manner, we introduce the expansions 

v ; _ v„ + s 2 v'i +... 



As - -x*+s As i + s 2 As 2 + ... 

We will write 

Ax k -Au k + As k for£:>l. (31) 

Expand the anharmonic term as 

f(X + sAxi +s 2 Ax 2 + ...) = s ft +s 2 f 2 + ..., (32) 

where terms in the r.h.s. depend on the Axj. Since / is as- 
sumed to have an overall order e or higher, the calculation 
of the term f k requires only the knowledge of Axj for j < k. 
Equations d24t , (|26| i and d3Tb then yield the recurrence rela- 
tions 



Au k (t) = 



1 



A u - Ag 
As k (t) = -Au k (0)e A >' - 



S[A u ,f k -t], 
1 



A u — A 



■S[As,f k ;t], 



(33) 



(34) 



Am = 



s Ami + s At/2 + . . . 



Ax k (t) = Au k {t) + As k (t), 
from which it can be finally obtained 

V* = (Au - A s )Au k (0). 
The recursion relations ([33} can be successively evaluated as 
written for k = 1,2,... up to any desired order. 

For example, for the anharmonic force corresponding to the 
generic one-dimensional potential dTOb with only cubic and 
quartic terms, expansion (l32l gives 

fx = -C3X 2 , 

fi — — C4X 3 — 2ct,X Ax\. 

It is then obtained 



Au\(f) 
Asi(f) 
Axi(t) 

Au 2 (t) 



C3 



A u - A s 

C3 



S[A u ,X 2 ;t], 



(S[Au,X 2 ;0]e^ + S[A s ,X 2 ;t]), 
A u ~ A s v ' 

— ^— (S [Au,X 2 - 0]e A >> -S[A u ,X 2 ;t] +S [A s , X 2 ; t]) .. 

A u - As v ' 

1 



A u - 
C'4 



S[Au,2ciXAxi +c A X i ;t\ 



Au - A s 
From (f34-b we have that 



S[A u ,X";t] 



2c\ 



(Au - A s ) 



2 S T [a u , X(t) (e A < T S [A u , X 2 ; 0] - S [A u , X 2 ; r] + S [A S ,X 2 ; t]) ; t] . 



V* = -c 3 S[A n ,X 2 ;0l 



in agreement with Eq. (|30l , and 



2c 



V* = - c 4 S [A U ,X 3 - 0] - -r-^yS r [Au, X(f) (e A * T S [A*, X 2 ; 0]-S [A u , X 2 ; r] + S [A s , X 2 ; t]) ; 0] . 



(35) 



(36) 
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FIG. 2. Critical velocity for a realization of the noise for a one- 
dimensional barrier with cubic anharmonicity, c 3 , for a> b = 1, y = 2, 
k B T= 1: 

Numerical simulation results (red crosses), harmonic ap- 
proximation ( 123) (gray horizontal line), perturbative results 
to first-order (123 1 -t- <T35T> (green straight line) and second- 
order P5)+p5ll+p6t (blue line). 



not only the position, but also the shape of the invariant man- 
ifolds, are stochastically time dependent. 

We also calculated the critical velocity numerically for a 
given realization of the noise. To this end, an ensemble of 
trajectories starting on the DS was propagated numerically. 
By recording which trajectories were reactive and which were 
not, the value of the critical velocity could be bracketed with 
high accuracy. For one fixed realization and for a potential 
with only a cubic anharmonic term, the perturbative expan- 
sion is compared to numerical results in Fig. [2] There is good 
agreement between perturbative and numerical results. Simi- 
lar figures are obtained for other realizations of the noise, thus 
leading to the same conclusion. Obviously, the size of the first 
and second order corrections, as well as that of the higher or- 
der terms that are omitted, varies among different realizations. 

In the special case that the anharmonic potential contains 
only a quartic term, the perturbation expansion results as an 
expansion in powers of e 2 , with the odd orders terms null. For 
the first two non-zero corrections, a similar calculation shows 
that 



Not surprisingly, the corrections (1351 and (1361 depend, 
through the function X, on the realization of the noise. This 
dependence reflects the fact that on an anharmonic potential 

I 



V* = -c 4 S[A u ,X 3 ;0], 



which is again consistent with Eq. < f30b . and 



(37) 



A u — A 



■S T [a u ,X 2 (t) (e' kT S [A a ,X 3 ; 0] - S [A a , X 3 ; r] + S [A S ,X 3 ; r]) ; 0] . 



(38) 



A comparison of the perturbative corrections (l37l and 
with numerical results is shown in Fig. [3] Again, this compar- 
ison confirms the accuracy of the perturbative results. 

The function X introduced in Eq. d29l plays a special role 
in the perturbation expansion because it represents the un- 
perturbed trajectory. To obtain a different perspective of this 
function, note that the critical velocity should depend only on 
the behavior of the stochastic force £ a (i) for t > 0, but not on 
the driving at earlier times: Once the initial conditions of a 
trajectory at t — are given, its future fate can only depend on 
the future noise. The separatrix between reactive and nonre- 
active trajectories must therefore also be determined by only 
the future noise. Yet the perturbation term in (Oil depends, 
via x*(f), on s*(f), which is given by past noise. 

If we split up the integration range of the S functional, we 
find that for t > 



s*(f) 



^V(0)+ f ' e' u '- T) 
Jo 



£*(t) dr. 



The integral in this expression depends only on noise for t > 0. 
The term including s*(0) contains all the dependence on the 
past, and it drops out when we form X{f). The variable X is 
the simplest modification of x* in which the dependence on 
the past has been removed. 



V. CORRECTIONS TO THE REACTION RATES 

A. General rate expressions 

In a one-dimensional model, the characteristic function Xt 
can be expressed in terms of the critical velocity as 



v x < V*. 



(39) 



In contrast to the TST approximation ((5), and in spite of its 
simplicity, the expression ( 1391 is exact. It allows to evaluate 
the average over initial conditions in Eq. © — the factor p ± in 
Eq. (0 being absent in one dimension — to find 



k — ( exp - 



2kuT 



(40) 



where only the average over the noise remains. This expres- 
sion was derived in Ref. for a harmonic barrier. It is now 
clear that the same expression holds also for anharmonic po- 
tentials if the critical velocity is suitably modified. Re- 
markably, no anharmonic corrections arise in the rate expres- 
sion d40b . 
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FIG. 3. Critical velocity for one realization of the noise for a 
one-dimensional barrier with quartic anharmonicity, C4, for cj b = 1, 
y = 2.5, k B T = 1: 

Numerical simulation results (red crosses), harmonic ap- 
proximation ( 123) (gray horizontal line), perturbative results 
to first-order (123l +d37ll (green straight line) and second- 
order @3)+@7J+@B) (blue line). 



If we have a perturbative expansion 



V* = V* + sV\ + e l W\ + 



(41) 



we can substitute into d40b and expand the exponential to ob- 
tain a series of rate corrections 



where 



K = Kq + EK\ + E K2 + ■ ■ ■ , 

K0 = (P)a , 

1 -t( pv *\ 



2k B T 

with the abbreviation 

n \ 



P = exp 



V, 







2k R T 



— exp 



(A u -A s ) 2 uV(0) 
2kuT 



(42) 

(43a) 
(43b) 

(43c) 
(44) 



We will now address the problem of evaluating the noise av- 
erages in Eq. ( |43l . 



B. Distorted correlation functions 

The corrections to the critical velocity that appear in the av- 
erages ( l43l are expressed in terms of the function X(t), which 
is in turn given in terms of the components u l {t) and s*(t) of 
the TS trajectory. They are Gaussian random variables whose 
correlation functions were evaluated in Ref. HH In the current 
notation and with 



k B Ty 



(45) 



2 A,t 

a e 1 , 



they read, for t > 0, as 

(sHt)sHo)) a 

(u t (t)u t (0)) a 

( u Ht)sHo)} a 



^crV' 1 "', 



= 0, 



(46a) 
(46b) 
(46c) 
(46d) 



To evaluate the corrections (|43l to the reaction rate, we 
need to calculate noise averages of the form (P(. . . )) a , where 
(...) indicates some expression in the functions «*(?) and 
s*(t). We will therefore assume that the expression (...) can 
be written as a function of finitely many random variables 
Z = (zi, . . .,z„) that follow a multidimensional Gaussian dis- 
tribution with zero mean and covariance matrix E, i.e., the 
matrix elements of S are cr,j = (ziz/j ■ As the first component 

we include the variable zi = m*(0), which plays a special role 
because it occurs in the factor P in Eq. d44b . 

Using d23l and setting p = (A u - A s ) 2 /k B T, we can write 



(P(...))a = 



1 



V(27r)"detE 
1 

V(2tt)" detE 



ze Z ^ zl2 e-P^\...) 



(...) 



/jgtfr 1 _ r tVt/2 

V detE V (27r)»det£ J K >' 



detZn 



(..•>o, 



detS 

where we have introduced the matrix 



(47) 



J = 



A 





and we have used (...) n to denote an average over a multidi- 
mensional Gaussian distribution with the modified covariance 
matrix En given by 

E ' = 2T 1 +pJ. 



From the observation 



1J 



0-11 

021 
<T„l 







we obtain (E7) = o~\{LJ. It is then easy to check that 

1 1 I/EWE -1 +pj) = /, 

\ 1+po-n ' 

the identity matrix. Therefore 



1 + pern 



■I/E 
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= Z + — pZ/Z, 



(48) 



where in the last step we have used the value given in Eq. (l46l l 
for o-n = (w* 2 (0)} . 

Furthermore, 



ZqZ- 1 =/- 



1 +po"i 



-Z/ 



is a lower triangular matrix whose diagonal elements, except 
for the (1,1) element, are all equal to 1. This observation 
makes it easy to evaluate 



det Z 
detZ 



det / 



P 



1 



1 + pern 
po"n 



-Z/ 



1 + pern 
_ ^ 
A s cJf 



where Eq. d46l i has again be used. 



(49) 



Substituting Eq. d49l in Eq. (l47b . we finally find 



(P(...)) a = — (...) . 



(50) 



For the components of the modified covariance matrix (PfFt we 
find 

which allows to obtain the moments of the distorted Gaussian 
distribution once the moments of the original Gaussian are 
known. In particular, \ZiZjJ Q = \ZiZjJ if either zi or Zj are 
uncorrelated with M*(0). 

Once the second moments of the distorted Gaussian distri- 
bution, i.e., the matrix elements of Zo, are known, Isserlis' 
theorem^! can be used to express higher-order moments in 
terms of second moments, e.g. 

<ZlZ2£3Z4>0 

= <ZlZ2>0 feZ4>0 + <ZlZ3>() feZ4>0 + <ZlZ4>0 (ZlZi\) ■ 

This expression contains a sum over all possible pairings of 
the four factors. Other even-order moments can be evaluated 
in a similar way, and the odd-order moments are zero. In this 
way, the modified averages of arbitrary polynomials can be 
calculated. 

The moments that will be required in the rate calculation 
can be obtained from these results; they are 



(«*(0)X(/)) 
{X(t)X(t')) 
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= (1-A> 



A u 



(52a) 



|l _ 2y 6 s + iij e~ 1 ^ +(l-fia) [e- 1 ^ + e^^') , (52b) 



with 



2A U 



2i s 



A u +A s 



C. Results for the one-dimensional potential 

With the help of Eq. (l50l l the leading term in the transmis- 
sion factor ( 143 at can be evaluated, giving 



An 

a> b 



(53) 



This is the famous Kramers result for the transmission 
factor. 17 

The perturbation expansion is set up in such a way that ef- 
fectively the noise carries a factor of s. The critical velocity 



Vj[ is linear in the noise. If Vj is one order e higher, it must be 
quadratic in the noise, and v\ cubic. Consequently, 



Kl = - 



1 A, 



Ml). 







is a third-order moment of the noise and must vanish. Simi- 
larly, all odd-order corrections to the transmission factor must 
be zero. According to the fluctuation-dissipation theorem (0, 
the noise carries a factor y/k^T, so that a perturbative expan- 
sion in powers of s corresponds to an expansion in powers of 
y/kftT. By contrast, Eq. (l42l is an expansion of the transmis- 
sion factor in powers of k%T because it has only even-order 
terms. 

The simplest rate correction can therefore be obtained from 
a quartic perturbation in the potential. We set C3 = 0, which 
makes Vf = 0, and calculate the rate correction that is linear 
in C4. Substituting Eqns. ( 1231 and ( 137] ) into ( 143bl i. it is found 
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that 



k 4 = 

K 2 



c 4 (A u - A s ) 
k B T 



S T [A u ,{PuH0)X\T)) a ,0]. (54) 



The average over the noise can be brought inside the S func- 
tional because the latter is shorthand notation for an integral. 
The remaining moment can be evaluated as 



(p^o)x\T)) a = ^-(uHo)x\T)) o 

= 3^( M *(0)X(r)) () (x 2 (r)) . (55) 



The modified correlation functions that are required here are 
given in Eq. (1521 , Equation (l55l l can thus be rewritten as a sum 
of exponentially decaying terms, for which the S functional in 
Eq. ( 1541 is easy to evaluate. This procedure yields 



c 3c 4 cr 4 (A u - A,) 2 3 
K 2 = TT~^ ; = 4 c 4 k B T 



4-k B Taj b A u 



(o>A s (A u - A s ) 2 



(56) 



This result agrees with the perturbative correction given in 
Ref s . [26M281 It can be rewritten as 



3 c 4 k B T (I - n 



4 w 4 \ l + ^2 



(57) 



in terms of the dimensionless parameter fi = Ko — /i u / w b that 
was used in Ref. [271 A comparison of Eq. d56l l with numer- 
ical results is shown in Figure |4] They confirm once more 
that the perturbative result is correct. The figure also shows 
the second-order correction in c 4 , which can be obtained in a 
similar way from Eq. (l38l . It reads 



3 

"32 



' c 4 k B T 



105// 8 + 830pt 6 + 1648/ + 770/x 2 + 87 
(1 -/)(3/ + 10/ + 3) ' 



(58) 



In the numerical example the second-order contribution is 
small, but Fig. @tb) shows clearly that the second-order per- 
turbative result is in better agreement with the numerical data 
than the first-order result. 

For a generic anharmonic potential that has a third-order 
term, the leading rate correction is quadratic in C3 and can be 
obtained from Eq. (I43cl > with the help of Eqns. d35l l and d36b . 
It reads 



\c\k B T (l-fi 2 \ 2 10/ +41^ 2 + 10 



K0 



1 +/ 



2fi 4 + 5ju 2 + 2 



(59) 



A comparison between Eq. d59t and numerical data is shown 
in Fig. Again, the agreement is excellent. 

If both cubic and quartic perturbations are present in the 
potential, then the second order contribution to the Kramers' 
transmission factor equals to the sum of expressions (IBTj 
and i 
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FIG. 4. Transmission factor, k, for a one-dimensional potential with 
quartic anharmonicity, c 4 , for a)^ = 3, k^T = 1. 

(a) a- as a function of the coupling strength c 4 for a value of the damp- 
ing y = 7. 

(b) Difference between k and its Kramers approximation, kq, as a 
function of y for c 4 =2: 

Numerical simulation results (red points), harmonic (Kramers) ap- 
proximation d53t (gray horizontal line), perturbative results to first- 
order, obtained from (I53t +(l56t (green line), and second-order ob- 
tained from (ED+{56}+OD (blue line). 




-0.3 -0.2 -0.1 0.1 0.2 0.3 



FIG. 5. Transmission factor for a one-dimensional potential with 
cubic anharmonicity, c 3 , with w b = 1, y = 2, k 3 T = 1: 
Numerical simulation results (red points), harmonic (Kramers) ap- 
proximation d53| l (gray horizontal line), perturbative results to 
second-order, obtained from (I53i +(f59l> (blue line). Notice that in 
this case the first-order correction is zero. 
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VI. THE TWO-DIMENSIONAL CASE 

So far, our discussion of the stochastic stable and unsta- 
ble manifolds and their use has been restricted to a one- 
dimensional model. Most problems of physical interest, how- 
ever, have several degrees of freedom. It is therefore crucial 
to show how the results obtained before can be generalized to 
higher dimension. We will carry out the generalization to two 
dimensions, which requires some extensions of the previous 
discussion. It will then be obvious that these techniques can 
equally be applied to systems in arbitrary dimension. 

We study a two-dimensional model whose dynamics is de- 
scribed by the Langevin equation ©. We denote the config- 
uration space coordinates as q — (x,y) and the correspond- 
ing velocities as q — (v v , v y ). The friction matrix F = yl 2 is 
assumed to be a scalar multiple of the 2x2 identity matrix, 
I 2 . By the fluctuation-dissipation theorem (0, this assumption 
implies that the x and y components of the fluctuating force are 
statistically uncorrected. For demonstration purposes we will 
use the anharmonic model potential 

U(x, y) = - l -u>l x 2 + io^y z + c x 2 y 2 (60) 

that has already been used in Refs.[32|and[33l The anharmonic 
perturbation in (1601 is of fourth order. In the terminology of 
the previous sections, the coupling parameter c is therefore of 
order e 2 , and rate corrections at first order in c are expected. 



A. Invariant manifolds in higher dimension 

In a two-dimensional setting, the phase space of the 
Langevin equation ^ is four-dimensional. It can be described 
with coordinates (x,y,v x ,v y ). As before, the harmonic ap- 
proximation of the dynamics around the barrier can be di- 
agonalized by introducing the coordinates u and s given in 
Eq. ( fT2b and coordinates z\ and z 2 defined by 



v y - A 2 y 



Z2 



y y -My 

Ai — A] 



(61) 



A\ - Az 

with the inverse transformation 

y = z\ + zi, v y = Am + A 2 zi- (62) 

The two additional eigenvalues 



lj2 = --( r± ^/y 2 - 4w 2 ) 



(63) 



are either real and negative or form a pair of complex conju- 
gates with negative real parts. 

The fluctuating force has two independent components 
%x,a(t) and £y, a (f), which determine the four components of 
the TS trajectory 



u\t) = 



1 



S[A u ,^ a ;t], 



s*(t) = ———S[A t ,£ ]Vt ;i], 

A a - A s 



zUt) = 



1 



S[A u ty, a ;tl 



A\ - A 2 

z\\t) = -—^—S[A2^ y , a -,t\ 
Ai-A z 



(64) 



that serves as a time-dependent coordinate origin. In the rela- 
tive coordinates 



Au = u — «*, As — s — s , 
Az\ - zi - z\ , Az 2 = Z2 - z\ 
the Langevin equation is written as 

fx(x,y) 



(65) 



Am = A u Au + 
As = A s As - 
Azi = /liAzi + 
Az 2 = A 2 Az2 - 



fx(x,y) 
A u - A s 
f y (x,y) 
A\ - A 2 
f y (x,y) 
A] — A? 



(66) 



where f x and f y denote the anharmonic parts of the mean 
force: 

~~fa = u b x + f x (x,y), 



= -(o y y + f y (x,y). 



_dU_ 

~dy 

The differential equations ( |6*6*] l are coupled by the conditions 

x — x" + Au + As, 

y = y* + Az\ + Az 2 . 

As in the one-dimensional case, the equations of mo- 
tion (l66*l l decouple and become time-independent in the har- 
monic limit, f x = f y = 0, and the relevant phase space struc- 
tures can easily be described in this case. Among the eigen- 
values in Eq. (l66l l. A u is positive, while the other three have 
negative real parts. Consequently, the TS trajectory has a one- 
dimensional unstable manifold and a three-dimensional stable 
manifold. The stable manifold separates reactive from non- 
reactive regions of phase space. The dimension of the unstable 
manifold, by contrast, is too low to separate distinct regions 
in the four-dimensional phase space. The invariant manifolds 
cannot therefore be used to distinguish trajectories with dif- 
ferent behaviors in the remote past, but the stable manifold 
can be used to predict the fate of a trajectory in the future. 
Thus, in arbitrary dimension the invariant manifolds provide 
precisely the diagnostic capabilities that are needed for rate 
calculations. 

We are particularly interested in trajectories that start on 
the DS x — 0. This is a three-dimensional surface with co- 
ordinates (v x ,y, v y ), embedded in the four-dimensional phase 
space. It intersects the three-dimensional stable manifold in 
a two-dimensional surface that separates reactive from non- 
reactive trajectories within the DS. We will call that two- 
dimensional surface the separatrix, and it depends on the real- 
ization of the noise. 
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FIG. 6. Schematic representation of the separatrix within the divid- 
ing surface x = 0. For a harmonic barrier the separatrix is a plane 
(gray in both panels), (a) For a weakly anharmonic barrier the sepa- 
ratrix can be parameterized by a function V*(y, v v ). Trajectories with 
initial condition v, > V*(y, v y ) are reactive, (b) If anharmonicities are 
strong, the separatix cannot be described by a single critical velocity, 
V*. 



On physical grounds, we expect a trajectory to be reactive 
if its initial velocity v x is sufficiently high. The critical ve- 
locity that separates reactive from non-reactive trajectories 
depends, in general, on the transverse coordinates y and v v . In 
the harmonic limit, the critical velocity is given by (l23l and 
is independent of these transverse coordinates. The separa- 
trix v x = V 1 is therefore a plane within the DS that is parallel 
to the y-Vy plane. When anharmonicities are taken into ac- 
count, the separatrix is deformed from this plane in a stochas- 
tically time-dependent way, as indicated schematically in Fig- 
ure |6ja). Nevertheless, we will still be able to describe the 
separatrix by specifying a critical velocity that depends on the 
transverse coordinates. In Section IVIBI a perturbative expan- 
sion for the function V*(y, v v ) will be developed. 

It is instructive to study the actual shape of the separatrix in 
a representative example. Figure |7]shows the critical velocity 
as a function of transverse coordinates for one realization of 
the noise for the two-dimensional model potential ( |60b . The 
critical velocity takes a maximum that is noticeably displaced 




FIG. 7. Critical velocity as a function of the transverse coordinates 
for one realization of the noise for the two-dimensional model poten- 
tial |60]l for <±> A . = 1, cj y = 1.5, y = 2, c = 0.2, k B T = 1. Contour 
spacing is 0.2 and the central contour value is -3.2. 

from the origin y = v y = 0. At the maximum, the critical 
velocity is closest to its harmonic value, which in this case is 
approximately -3.01. For all values of the transverse coordi- 
nates, the critical velocity is below the harmonic approxima- 
tion value. Moreover, it decays steeply away from the maxi- 
mum, so that deviations from the harmonic approximation are 
large for most values of the coordinates. As the critical veloc- 
ity appears in the exponent in the rate formula (|40l — which 
will be generalized to higher dimension in Eq. ( |73l — , it is 
expected that anharmonic effects on the critical velocity leads 
to large rate corrections. 

If the barrier is strongly anharmonic it cannot be guaran- 
teed, in general, that the separatrix can be parameterized by 
the transverse coordinates y and v y . In a situation as that indi- 
cated in Fig.|6tb), the separatrix is described by a multivalued 
function of the transverse coordinates. It cannot be charac- 
terized by a single critical velocity. As expected. trajectories 
at low v x are nonreactive, and those at somewhat larger v x 
are reactive. However, at certain values of y and v v , there is 
an interval at yet higher v x that also contains nonreactive tra- 
jectories. A scenario like this obviously requires very strong 
anharmonic effects, and this is only be achieved for large val- 
ues of the transverse coordinates. But at these conditions, it 
is doubtful whether a TST-like treatment with a single rate- 
determining saddle point is appropriate at all. We will there- 
fore neglect this possibility and assume the existence of a sin- 
gle critical velocity. 



B. Determination of the stable manifold 

As a basis for the perturbative expansion, we formally solve 
the differential equations d66| i in terms of S functionals by 

Au(t) = — !— S[A u ,f x (x,y);t], 

As(0 = As(0)e As ' - S [As,f x (x,y); t], 
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Azi(t) = Azi(0)e 



Az 2 (f) = Az 2 (0K 



1 



A\ - A 2 
1 

A t - A2 



s[A u f y (x, y y,ti 

S[A 2 ,Mx,y);t]. (67) 



These integral equations are entirely analogous to Eqns. (l24l 
and d26l l. and they are coupled by 

x = x* + Au + As, 
y = y + + Azi + Az 2 . 

A trajectory satisfying (l67l i automatically lies on the stable 
manifold. To find the critical velocity, Eqns. (|67T > needs to be 
solved under the condition that the trajectory starts in the DS 
x — and at the prescribed transverse coordinates y(Q) and 
v,(0). 

We will solve Eqns. ( f6Tb by an iterative procedure as 
in (I33t . As before, the initial condition As(0) must be adapted 
in every step in order to enforce the condition x(0) = 0. By 
contrast, the transverse initial conditions Azi (0) and Az 2 (0) are 
fixed once and for all by imposing the condition that 

y(0)=y*(0) + Azi(0) + A?2(0) t 
v v (0) = v*(0) + AiAzi(0) + A 2 Azi(0) 

take the desired values. The critical velocity is finally obtained 
from Eq. (127) . 

Our perturbation expansion is centered around the har- 
monic approximation to a trajectory on the stable manifold, 
given by Eq. (f29]) 



X(t) = xHt) - x*(0)e 



and 



Y(t) = y*(f) + Azi (O)e- 11 ' + Az 2 (0)e i2f . (68) 
The latter can be split according to 

Y(t) = Y a (t) + Y x (t) (69) 

into one part 

Y a (t) = yHt) - z\(Q)e Alt - z\{Q)e^' 

that depends on the realization of the noise but not on the ini- 
tial conditions, and another 



Y x (t) = zi(0)e M ' + z 2 (0)e 

that depends on the initial conditions but not on the noise. 

We will now apply the general theory to the model poten- 
tial (160b . Our aim is to expand the coordinates 



y{t) = Y{t) + cAy l {t) + c 1 Ay 1 {t) + ... 

in powers of the anharmonicity parameter c. For expansions 
of other quantities, such as 

V* = V* + cV* + c 2 y| + . . . , 
a similar notation will be used. The anharmonic forces are 
given by 

f x = -2c xy 1 

= -2c XY 2 - 2c 2 (Y 2 A Xl + 2XY Ay, ) + . . . , 
fy = ~2cx 2 y 

= -2cX 2 Y - 2c 2 (2XYA Xl + X 2 A yi ) + .... 

In the first step of the iteration we find 
1 



Aui(t) 



A u — A. 



S[A u J x y,i\ 



■S[A u ,XY 2 ;t], 



(70) 



where f x n is the coefficient of f(x) of order c". From Eq. (l70l 
we get 



V* = (A a - A s )Au l (0) 
= -2S[A U ,XY 2 ;0]. 



(71) 



The remaining coordinates need only be calculated if the 
second-order correction for the critical velocity is desired. We 
then obtain 

A Sl (0 = -A Ml (0y»' + — ?— S[A s ,XY 2 ;t], 
A u - A s 

Azi(t) = —-^—S[A u X 2 Y;t], 
A\ - Az 

Az 2 (t) = +—?—S[A 2 ,X 2 Y;t]. 
A\ - A 2 

Finally, with the aid of 

Ax\ — Au\ + Asi, 
Ay\ = Azi + Az 2 , 

we can calculate 



Au 2 (t) 



S[A u ,f x x,t]. 



x(t) = X(t) + c Axi(t) + c Ax 2 (t) + . 



The resulting expression reduces to 
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V| = -45 T 



Jk, y-zj(s [A u ,XY 2 ;0]e A > T - S [a u ,XY 2 ;t]+ S [a s ,XY 2 ;t]} 



+ 2 



X(t)Y(t) 
A\ - A<i 



(s [a 2 ,X 2 Y;t] -S [a u X 2 Y;t]\, 



(72) 
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FIG. 8. Critical velocity for one realization of the noise for the 
two-dimensional model potential ( 1601 ) with co x = 1, co y = 1.5, y = 2, 
k%T = 1, for an initial condition y = 0, v y = 0. 
Numerical simulation results (red crosses), harmonic ap- 
proximation d23t (gray horizontal line), perturbative results 
to first-order J23t +<f7Tb (green straight line) and second- 
order d23}+([7B+<lZ2ll (blue line). 



Figure [8] shows the value of the critical velocity for one real- 
ization of the noise for the two-dimensional model potential 
( f60b as a function of the coupling strengh, c, for the initial 
condition y = 0, v y = 0. It is compared to perturbative re- 
sults up to second order. As it can be seen, our perturbative 
results agree very well with those obtained numerically, thus 
showing the efficiency of our method. To further analyze the 
performance of our method, we show in Fig. |9]the difference 
between the numerically calculated critical velocity and the 
value obtained with our perturbative expansions for different 
values of the transverse coordinates, where it is clearly seen 
that it sensibly reduces as the order of the perturbation is in- 
creased. 



C. Reaction rate expressions 

The simple expression (l40l > for the transmission coefficient 
in terms of the critical velocity can easily be generalized to 
higher dimension. To achieve this, we start again from Eq. (O. 
Note first that in the denominator of Eq. (|7]i the average over 
the transverse coordinates has no effect since the TST approxi- 
mation to the characteristic function does not depend on them. 
In the numerator, we use again the form ((39) of the character- 
istic function and carry out the average over v A as before, to 



exp 



y*2 
2AVT 



(73) 



In this expression the average over the transverse coordinates, 
which is indicated by subscript JL, cannot be carried out imme- 
diately because the critical velocity depends on the transverse 
coordinates. 

Equation ( 1731 represents the simplest conceivable general- 
ization of Eq. d40b . It is remarkable that no modifications, be- 
yond the additional average over the transverse coordinates, 
are required. This is only possible because no anharmonic 
corrections are required for the denominator in Eq. (0. 

In the case of the model potential d60b . the distribution (O 
of the transverse coordinates is given by 



P±(y, v y ) = ~ exp 



v: + 



2 2 \ 



2kuT 



(74) 



i.e., it is a Gaussian distribution. The functions X and Y will 
then both have a Gaussian distribution, which allows us to 
evaluate the rate corrections by the method of Sec. IV B| For 
any expression involving w*(0), X and Y, we write 



(P(...)) a± = —(...) 0± 



(75) 



as in Eq. (|5Q| ). The average over the initial conditions is not in- 
volved in the transition from the noise average to the distorted 
average with correlation function (IBTj , because the noise and 
the initial conditions are uncorrected. 

Once we have a perturbative expansion of the critical ve- 
locity of the form d4TT l. expressions (l43l can be used for the 
expansion of the transmission factor. The only required modi- 
fication being to replace noise averages by averages over noise 
and the transverse coordinates. 

Assuming a general anharmonic potential of the form 

C/(0 ) y) = ^y 2 + C/ anh (y), 

where U al ±(y) contains terms at least of third order in y, i.e. 
at least of first order in the expansion parameter e, it can be 
treated perturbatively in the current framework. The distri- 
bution function of the transverse coordinates can then be ex- 
panded as 



pAy,Vy) = | exp 



I 2 , 2 i \ 

vj + gjft~ ^ 
2k B T 
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FIG. 9. Difference between numerically calculated critical velocity 
and perturbative expansions. Noise sequence and parameter values 
as in Fig. [7] (a) Harmonic approximation, (b) First-order perturba- 
tion theory, (c) Second order perturbation theory. Contour spacing 
is 0.05 in (a), 0.005 in (b) and (c). Note that the color scale is also 
stretched by a factor 10 in (a). 

x (l + eai(y) + s 2 a 2 (y) + . . .) (76) 

with suitable coefficients a, that are polynomials in y of de- 
gree at most i. We assume that the partition function Z in 
Eq. d76l ) is the same as in the Gaussian distribution d74l i, and 
any corrections to the partition function that arise from the 
anharmonicity of the potential have been included in the ex- 
pansion coefficients a,{y). 

Using symbol 1L to denote an average over the Gaussian 



distribution (l74l i of initial conditions, we can write 




The expansion fiTT i of the critical velocity then allows us to 
expand the exponential, thus obtaining 

K — Kq + SK\ + e 2 K2 + . . . 

with 

Ko = (P)aA . (77a) 
Kl = 'k^f ( PV W«ii + <^WU , (77b) 

(77c) 

wehre again the abbreviation (l44t has been used. The remain- 
ing averages are Gaussian averages that can be evaluated, as 
before, by first converting the noise average into a distorted 
Gaussian average via d75l l, and then using Isserlis' theorem. 

Because the factor P is independent of the initial conditions, 
we obtain from (177 at 

*o = <n, = -> 

the Kramers result. Similarly, the expressions (Pcii(y)) a]L , that 
occur in all correction terms, can be simplified to 

<Pa,-Cy)> oJ1 = {P) a (ai(y)) A = — (a®))* . 

D. Correlation functions 

To evaluate corrections to the transmission factor in 
Eq. ( TTTl i using Isserlis' theorem, the correlation functions 
(wiW2)oil, where w\ and wo are one of u*(0), X(t), Y(t), and 
y(Q), are needed. (The initial condition y(Q) was written with- 
out its time argument in Sec. lVICl For the sake of clarity we 
will now include it again.) 

Because the x and y components of the fluctuating force are 
uncorrelated, all correlation functions involving one of either 
m^O) or X(t) and one of either Y(t) or y(Q) must vanish. Fur- 
thermore, since u*(0) and X(t) do not depend on initial condi- 
tions, 

(u*(0)X(t)) QjL = (uHO)X(t)) o 

and 

(X(t)X(t')) 0jL = (X(t)X(t')) 
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are given by Eq. 

Concerning the initial conditions, it can be read off from the 
distribution function ( l74t that 

<^U = f- <™> 

(The average over the distorted noise distribution does not 
have any effect.) We can also see that 

(v y (0) 2 ) 0]L = k B T and (y(0) v y (0)) QR = 0. (79) 

These results further yield 

(y(0) Y(t)) 0A = (y(0) Y x (t)) 0]L 

(A l e A i'-A 2 e A < 1 ). (80) 



<y(0)zi(0)W 
k B T 



Finally, the autocorrelation function of Y(t) can be decom- 
posed, with the help of the split d69D . into 

(Y(t)Y(t')) 0]L = (Y a {t)Y a {t')) a + (Y±(t)Y ± (t')) ]L (81) 
because 

(Y a (t)Y ± (t')) 0]L = (Y a (t)) (Y x (t')) ± = 

and 

(Y a (t)Ya(f)) = (Y a {t)Y a {t')) a . 

To evaluate the first term in Eq. (ISTT i. the correlation function 
of the components z?(f) of the TS trajectory, given in Ref.l3ll 
are needed. The second term can be evaluated with the help 
of Eqns. d78l and (1791 1. Finally, one arrives to the following 
simple result 

(YW(t')) 0iL = T ^- 2 + (82) 



u>\,-X\ 



o£ - xl 



With that we have found all correlation functions that we will 
need to calculate the rate corrections. 



E. Rate corrections 

Let us now derive an expansion of the transmission factor 
for the case of the anharmonic model potential i 



K — Kq + C K\ + C K2 + ■ ■ 



(83) 



in powers of the coupling parameter c. As discussed earlier, 
this corresponds to an expansion in powers of e 2 , and the rate 
formulas (|77| | with afy) = can be used. 
The first correction term is 



' Oil 



K, =-— —(V^) 



(84) 



The remaining average can be simplified to 

{ U H0)X(t)Y 2 (t)) Q]l = ( U HO)X(t)) o (y 2 (t)} 



Oil 



The results of Sec. IVIDI give a sum of exponentially decay- 
ing terms for this expression, so that the S functional can be 
evaluated as in the one-dimensional case. In terms of the di- 
mensionless parameters fi = kq = X u /ci>\„ that was already 
used above, and v = coy/cob the rate correction reads 



Kl = — 



jk B T /u 2 



col ( 1+ J") 2y2 ' 



(85) 



The second-order correction can be obtained in a similar way. 
After tedious calculations, one finally arrives at 



J 



K2 = 



n(k B T) 2 l 96(m 2 -1) 2 6 
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9(at - 1)(3a* 4 + 8// 2 + 1) 


6tol W 2 + 1) V - 4v 2 - 2) {fi 2 + l)(v 2 + 1) 


(2// 2 + l)(3// 2 +4v 2 +6) 


{H 2 + 1)3 V* 



96(// + 2/r-l) y u 6 



(/a 2 + l) 2 (ju 2 + 2)(// 4 - 2^ 2 (2v 2 + 1) - 8) in 2 + 1)V 4 - 2// 2 (2v 2 + 1) - 3) 

192/ 2(16/* 12 - 24// 10 - 139^ 8 - 75/i 6 + 11^ + 1 1 \yr + 34) \ 



(2// 6 + 7/ + 7^ 2 + 2)(^ 2 (4v 2 + 6) + 3) 



(yu 2 + 1) 4 (2// 4 + 5ju 2 + 2)v 2 



(86) 



A numerical example is shown in Fig.[l0] The second-order 
corrections to the transmission coefficient are small, so that a 
large number of trajectories needs to be included in the nu- 
merical calculation of the rate. Nevertheless, it is clear that 
the perturbative expressions (l85"t and (l86t describe the rate 
correctly. 



VII. CONCLUDING REMARKS 



TST and related schemes have been widely used for rate 
calculations for a long time. For reactions that occur in so- 
lution, recrossings of the DS pose a major difficulty in such 
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FIG. 10. Transmission factor for the two-dimensional model poten- 
tial J601 > as a function of coupling strength, c, for oj b = 1, u y = 0.5, 
k B T= l,r= 1: 

Numerical simulation results (red points), harmonic (Kramers) ap- 
proximation < |53t > (gray horizontal line), perturbative results to first- 
order, obtained from (153> +ll85l (green line), and second-order ob- 
tained from (|53}+([85j+([86j (blue line). 



calculations. Many approaches try to overcome this problem 
by choosing the DS judiciously. By contrast, the method de- 
veloped here is insensitive to the choice of this surface. The 
simplest choice of DS, which was taken here, also leads to the 
simplest calculation of the critical velocity. The use of a dif- 
ferent DS would require a redefinition of the critical velocity 
to describe its intersection with the stable manifold, but this 
can be achieved with only minor modifications to the iteration 
procedure for the critical velocity. After that, any DS that lies 
within the barrier region would give the same rate. 

This independence of the DS is achieved by two crucial 
features of our method. First, the dynamics are described in 
phase space, rather than in configuration space, and modern 
geometric methods are used in our study. Second, we focus 
on invariant geometric structures that are determined by the 
dynamics, rather than in structures, such as the DS, that are 
arbitrarily imposed by the researcher. The present results in- 
dicate that similar results apply to reactive systems that are 
coupled to their environments, i.e., TST should focus on in- 
variant structures in phase space. 

The focus of the present paper has been on analytic pertur- 
bation theory for the rate corrections on an anharmonic bar- 
rier. The different steps of this calculation have different levels 
of complexity. The critical velocity, which encodes the loca- 
tion of the invariant manifold, is very simple to calculate with 
the iteration scheme described here. Moreover, it can easily 
be extended to higher orders. By contrast, the evaluation of 
the averages that yield the rate corrections is laborious. While 
straightforward in principle, it requires the calculation of a 
large number of exponential integrals something that, even 
for some of the results presented here, is only feasible with 
the help of a computer algebra system, Mathematical in our 
case. 

The crucial step that sets the current method apart from ear- 
lier algorithms is the calculation of the stable manifold and the 
critical velocity. It is encouraging, therefore, that this most 



important step of the calculation is also the easiest. This ob- 
servation further suggests that to obtain an efficient algorithm 
to compute rates, the calculation of the stable manifold should 
be combined with numerical methods for the computation of 
averages. We will report on such combinations in a forthcom- 
ing publication. 
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